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ABSTRACT 

In this paper a finite element method involving Petrov-Galerkin method with quintic B-splines as basis functions 
and septic B-splines as weight functions has been developed to solve a general sixth order boundary value problem with a 
particular case of boundary conditions. The basis functions are redefined into a new set of basis functions which vanish on 
the boundary where the Dirichlet and the Neumann type of boundary conditions are prescribed. The weight functions are 
also redefined into a new set of weight functions which in number match with the number of redefined basis functions. The 
proposed method was applied to solve several examples of sixth order linear and nonlinear boundary value problems. The 
obtained numerical results were found to be in good agreement with the exact solutions available in the literature. 

KEYWORDS: Absolute Error, Petrov - Galerk in Method, Quintic B-Spline, Septic B-Spline, Sixth Order Boundary 
Value Problem 

INTRODUCTION 

In this paper, we consider a general sixth order linear boundary value problem 

a 0 (x) y {6) (x) + a x (x) y {5) (x) + a 2 (x) y (4) (x) + a 3 (x) y"\x) + a 4 (x) y\x) 
+ a 5 (x)y'(x) + a 6 (x)y(x)=b(x), c<x<d 

subject to boundary conditions 

y(c)=A 0 , y(d)=C 0 , y'(c) =A U y\d) = C,, y"{c) =A 2 , y'\d) =C 2 (2) 

where A 0 , C 0 , A h C\, A 2j C 2 are finite real constants and a 0 (x), a\(x), a 2 (x), a 3 (x), a^x), a s (x), a 6 (x) and b(x) are all 
continuous functions defined on the interval [c, d]. 

The sixth order boundary value problems occur in astrophysics [1]. Chandrasekhar [2] determined that when an 
infinite horizontal layer of fluid is heated from below and is under the action of rotation, instability sets in. When this 
instability is as ordinary convection, the ordinary differential equation is sixth order. The existence and uniqueness of the 
solution for these types of problems have been discussed in Agarwal [3]. Finding the analytical solutions of such type of 
boundary value problems in general is not possible. Over the years, many researchers have worked on sixth-order 
boundary value problems by using different methods for numerical solutions. Wazwaz [4] developed the solution of 
special type of sixth order boundary value problems by using the modified Adomian decomposition method and he 
provided the solution in the form of a rapidly convergent series. Huan [5] presented variational approach technique to solve 
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a special case of sixth order boundary value problems. Noor et al. [6] presented the variational iteration principle to solve a 
special case of sixth order boundary value problems after transforming the given differential equation into a system of 
integral equations. Ghazala and Siddiqi [7] presented the solution of a special case of sixth order boundary value problems 
by using non-polynomial spline functions. Siddiqi et al. [8], Siddiqi and Ghazala [9] developed quintic spline funtions and 
septic spline functions techniques to solve a special case of linear sixth order boundary value problems respectively. 
Lamnii et al. [10], kasi viswanadham and Showri raju [11] developed septic spline collocation and quintic B-spline 
collocation method are used to solve sixth order boundary value problems respectively. Loghmani and Ahmadinia [12] 
used sixth degree B-spline functions to construct an approximation solution for sixth order boundary value problems. 
Waleed [13] presented Adomian decomposition method with Green's function to solve a special case of sixth order 
boundary value problems. Kasi Viswanadham and Murali krishna [14] developed septic B-spline Collocation method to 
solve a special case of sixth order boundary value problems. Kasi Viswanadham and Sreenivasulu [15] developed quintic 
B-spline Galerkin method to solve a general sixth order boundary value problems. So far, sixth order boundary value 
problems have not been solved by using Petrov-Galerkin method with quintic B-splines as basis functions and septic 
B-splines as weight functions. This motivated us to solve a sixth order boundary value problem by Pertrov-Galerkin 
method with quintic B-splines as basis functions and sextic B-splines as weight functions. The solution to a nonlinear 
problem has been obtained as the limit of a sequence of solution of linear problems generated by the quasilinearization 
technique [16]. Finally, in the last section, the conclusions are presented. 

JUSTIFICATION FOR USING PETROV-GALERKIN METHOD 

In Finite Element Method (FEM) the approximate solution can be written as a linear combination of basis 
functions which constitute a basis for the approximation space under consideration. FEM involves variational methods like 
Rayleigh Ritz method, Galerkin method, Least Squares method, Petrov-Galerkin method and Collocation method etc. In 
Petrov-Galerkin method, the residual of approximation is made orthogonal to the weight functions. When we use 
Petrov-Galerkin method, a weak form of approximation solution for a given differential equation exists and is unique under 
appropriate conditions [17, 18] irrespective of properties of a given differential operator. Further, a weak solution also 
tends to a classical solution of given differential equation, provided sufficient attention is given to the boundary conditions 
[19]. That means the basis functions should vanish on the boundary where the Dirichlet type of boundary conditions are 
prescribed and also the number of weight functions should match with the number of basis functions. Hence in this paper 
we employed the use of Petrov-Galerkin method with quintic B-splines as basis functions and septic B-splines as weight 
functions to approximate the solution of sixth order boundary value problem. 

DESCRIPTION OF THE METHOD 
Definition of Quintic B-Splines and Septic B-Splines 

The quintic B-splines and septic B-splines are defined in [20-22]. The existence of quintic spline interpolate s(x) 
to a function in a closed interval [c, d] for spaced knots (need not be evenly spaced) of a partition c = x 0 < Xj <...< x„.j < 
x n = d is established by constructing it. The construction of s(x) is done with the help of the quintic B-splines. Introduce ten 
additional knots x_ 5 , x_ 4 , x_ 3 , x_ 2 , x. u x n+i , x n+2 , x n+3 , x n+4 and x n+5 in such a way that 

X_5<X.4<X. 3 <X.2<X. 1 <Xo and X n <X n+ i<X n+ 2<X n+3 <X n+ 4<X n+5 . 

Now the quintic B-splines B t {xys are defined by 
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BXx) = 



— — , x£l L-*,_ 3 , x i+il 

r=i-3 7T{X r ) 



0 , otherwise 
where (x r - x) 5 + 



(x-x) , if x r >x 
0, ifx r <x 



!+3 

and K{x) = Y\ ( x ~ x r ) 

r=i-3 

where {_B_ 2 (X), B.j(x), B 0 (x), B 1 (x),...,B„. 1 (x), B„(x), B n+] (x), B n+2 (x)} forms a basis for the space Sj(;z) of quintic 
polynomial splines. Schoenberg [22] has proved that quintic B-splines are the unique nonzero splines of smallest compact 
support with the knots at 

X_5<X_4<JC_3<X_2<JC_1 <^Xq<^X\ < K.X n .\ <^X n <^X n+ i<^X a+ 2 < ^X n+ ^<^X n+ ^<^X n+ ^. 

In a similar analogue septic B-splines /?,(jc)'s are defined by 



RAx) = 



Z {x r x) + 

r=i-4 KyX r ) 

0, otherwise 

7 [(* r -*)\ if x r ^ x 



where (x - 



[0, if x ,.^ x 



j+4 

and ^(x) = |~[ (x-X r ) 

r=!-4 

where {/?_/x,), R.2(x), R.i(x), Ro(x), Ri(x),,..,R„.i(x), R„(x), R„+i(x), R n+2 (x), R„+3(x)} forms a basis for the 

space 5 7 (^") of septic polynomial splines with the introduction of four more additional knots x_ 7 , X-e, x n+6 , x n+1 to the 

already existing knots x_ 5 to x n+5 . Schoenberg [22] has proved that septic B-splines are the unique nonzero splines of 
smallest compact support with the knots at 

X_7<X_(,<X_s<X_4<X_3<X_2<X_i<Xo<X\<. . . <X n .i<X n <X n+ i<X n+2 <X n+ ^<X n+ 4<X n+ ^<X n+ (,<X n+ - l . 

To solve the boundary value problem (1) subject to boundary conditions (2) by the Petrov-Galerkin method with 
quintic B-splines as basis functions and septic B-splines as weight functions, we define the approximation foiy(x) as 

n+2 

y(x)=^a j B j (x) (3) 

where a/s are the nodal parameters to be determined and B/x)'s are quintic B-spline basis functions. In 
Petrov-Galerkin method, the basis functions should vanish on the boundary where the Dirichlet type of boundary 
conditions are specified. In the set of quintic B-splines {B. 2 (x), B.j(x), B 0 (x), B ] (x),...,B„_ 1 (x), B n (x), B n+1 (x), B n+2 (x)}, the 
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basis functions B. 2 (x), B A (x), B 0 (x), B Y (x), B 2 (x), B n . 2 (x), B nA (x), BJx), B n+1 (x) and B n+2 (x) do not vanish at one of the 
boundary points. So, there is a necessity of redefining the basis functions into a new set of basis functions which vanish on 
the boundary where the Dirichlet type of boundary conditions are specified. When the chosen approximation satisfies the 
prescribed boundary conditions or most of the boundary conditions, it gives better approximation results. In view of this, 
the basis functions are redefined into a new set of basis functions which vanish on the boundary where the Dirichlet and 
Neumann type of boundary conditions are prescribed. The procedure for redefining of the basis functions is as follows. 

Using the definition of quintic B-splines, the Dirichlet and the Neumann boundary conditions of (2), we get the 
approximate solution at the boundary points as 



\ = y(c) = y(x Q ) = £ (ZjBjiXo) 



j=-2 



(4) 



n+2 



C 0 = y(d) = y(x n ) = X OC.B.{x n ) 



j=n-2 



(5) 



A = /(c) = y'(x 0 ) = £ ajB'j(x Q ) 



j=-2 



(6) 



n+2 



C l = y\d) = y\x n )= X OC.B'.{x n ) 

j=n-2 



(7) 



Eliminating a_ 2 , a_i, a n+ i and a n+2 from the equations (3) to (7), we get 



y(x)=w(x) + Y j a j Q j (x) 

J=0 



(8) 



where 



(9) 



w 1 (x) = -^-B_ 2 (x)+-^-B n+2 W 
B_ 2 (x Q ) B n+2 (x n ) 



(10) 



7 = 0,1,2 



j = 3,4,...,n-3 



P'(x ) 

Pj(x)- J " P n+l (x), j = n-2,n-\,n 



(ID 
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B.(x 0 ) 

-B_ 2 (x), y = -1,0,1,2 
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Pj(x)= 



B- 2 (x 0 ) 



Bj(x)- 

Bj(x), 

B.(x n ) 

B i^~~^ — — fi „ +2 W' j = n-2,n-l,n,n + l 



j =3,4,...,n— 3 



(12) 



The new set of basis functions in the approximation y(x) is {Q/x), j=0,l,...,n}. Here w(x) takes care of given set 
of Dirichlet and Neumann type boundary conditions and Q/x)'s and its first order derivatives vanish on the boundary. In 
Petrov-Galerkin method, the number of basis functions in the approximation should match with the number of weight 
functions. Here the number of basis functions in the approximation for y(x) defined in (8) is n+1 , where as the number of 
weight functions is n+1. So, there is a need to redefine the weight functions into a new set of weight functions which in 
number match with the number of basis functions. The procedure for redefining the weight functions is as follows: 

Let us write the approximation for v(x) as 



(13) 



7=-3 



where R/x)'s are septic B-splines and here we assume that above approximation v(x) satisfies corresponding 
homogeneous boundary conditions of the given boundary conditions (2). That means v(x) defined in (13) satisfies the 
conditions 



v(c) =0,v(d) = 0, v (c) = 0, v\d) = 0, v"(c) = 0, v\d) = 0 

Applying the boundary conditions (14) to (13), we get the approximate solution at the boundary points as 

v(c) = v(x 0 ) = f d j3 j R j (x 0 ) = 0 

7=-3 

v(d) = v(x n )= ]T j3 j R J (x n ) = 0 

j=n-3 

v'(c) = v\x 0 )=f j /3 j R' j (x Q ) = 0 

7=-3 

v\d) = v\x n )= X J3jRj(x n ) = 0 

j= n -3 

v"(c)=v"d„)=x^;u)=o 



(14) 



7=-3 



n+3 



\d) = v'\x n )= XM"^) = 0 

j=n-3 



(15) 



(16) 



(17) 



(18) 



(19) 



(20) 
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Eliminating p_ 3 , P_ 2 , p. i, P„+i, P„+2 and p n+3 from the equations (13) and (15) to (20), we get the approximation 
for v(x) as 



v(*) = £#V,(*) 
where 



Vj(x) = 



r.w-^fr^x), ; = o,i,2,3 



T -i( x o) 



Tj(x), 



7=4,5,..., n-4 



T j (x) = 



T'\x ) 

Tj (x) - ' " T n+1 (x), j = n-3,n-2,n-l,n 

s M)-7^r^ s -2(x)> ; = -i,o,i,2,3 



Sj(x), 



7 = 4,5,...,n-4 



^■W- ^ " s S n+2 (x), j = n-3,n-2,n-l,n,n + l 



(21) 



(22) 



(23) 



Sj(x) = 



R ( x ) 

R j W - TTT^T W • J = -2,-1,0, 1, 2, 3 



Rj(x), 



7=4,5,...,n-4 



7? ,. (x) i? . , (x) , j = n — 3, n - 2, n - 1 , n, n + 1, n + 2 



(24) 



Now the new set of weight functions for the approximation v(x) is {V/x), j=0,l,...,n}. Here V/x)'s and its first 
and second order derivatives vanish on the boundary. 

Applying the Petrov-Galerkin method to (1) with the new set of basis functions {Qj(x), j=0, 1,..., n} and the new 
set of weight functions {V/x), 7=0,1, ...,n}, we get 



{[^(x)y^x)+^(x)y 5) (x)+a 2 (x)y^x)+a 3 (x)/ / (x)+ fl4 (x)Ax)+ fl5 (x)/W 
+a (i (x)y(x)]V i (x) dx= j b{x)V i (x)dx for i=0, l,...,a 



(25) 



Integrating by parts the first two terms on the left hand side of (25) and after applying the boundary conditions 
prescribed in (2), we get 



\a»{x)V i {x)y«\x)dx = -^- J [a 0 {x)V i {x)\ Q+^rkWTOl ^+ \^ J [a,{x)V i {x)\y"{x)dx (26) 



d' r , . X C d 4 



dx 



dx" 
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x " x " d 3 

f a, {x)V i (x) y {5) (x)dx = - [ — [a, (x)V i (jc)] y"(x)dx (27) 

Substituting (26) and (27) in (25) and using the approximation for y(x) given in (8), and after rearranging the 
terms for resulting equations, we get a system of equations in the matrix form as 

Aa = B (28) 
where 

A = [a ij ]; 

a. . = P {a 2 (x)V ( (x)fi/ 4) U) + a,(x)V i .(x)Q"\x) + [-^-[a 0 (x)V;(*)] - ^[ fll (x)V,(x)] 

Jx » ax ax 

+a 4 (x)V i (x)]Qj(x) + a 5 (x)V, (x)Q' j (x) + a 6 (jc)V, . (x) }dx 

for i=0, 1, ... ,n; j=0, 1, ... ,n. (29) 
B = [£,]; 

b, = P {fc(*)V; (jc) - {a 2 (x)V i (x) w (4) (x) + a 3 (x)V i (x) w'"(x) + K (*W (*)] ~\(h (W to] 

<2X <3X 

^3 ^3 

+a 4 (x)V;. (x)]w"(x) + a 5 (x)VJ (x) w'{x) + a 6 (x)V; (x)w(x) } }dx + — [a 0 (x)V; (x)] C 2 — — [a 0 (x)V t (x)l A, 

ax dx " 

for i=0, 1, ...,«. (30) 

and a = [a 0 a l ...a n f . 

PROCEDURE TO FIND THE SOLUTION FOR NODAL PARAMETERS 

A typical integral element in the matrix A is 

n-i 

«7=0 

where / = f ' ' v (x)r (x)Z(x)dx an d fj{x) are the quintic B-spline basis functions or their derivatives. V ; (x) 

are the septic B-spline weight functions or their derivatives. It may be noted that /„, = 0 if 
(x i ._ 4 ,x. +4 )n(x._ 3 ,x. +3 )n(x m ,x m+1 ) = 0- To evaluate each/ m , we employed 7-point Gauss-Legendre quadrature formula. 

Thus the stiffness matrix A is a thirteen diagonal band matrix. The nodal parameter vector CC has been obtained from the 
system A(X = H using the band matrix solution package. We have used the FORTRAN -90 program to solve the 
boundary value problems (1) - (2) by the proposed method. 

NUMERICAL RESULTS 

To demonstrate the applicability of the proposed method for solving the sixth order boundary value problems of 
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the type (1) and (2), we considered three linear and three nonlinear boundary value problems. The obtained numerical 
results for each problem are presented in tabular forms and compared with the exact solutions available in the literature. 

Example 1: Consider the linear boundary value problem 

_4y 4 > +2y" + xy = (5 + 2x-x 2 )e x , 0<x<l (31) 

subject to 

y(0) = 1, y(l) = 0, /(0) = 0, y\\) = -e, y"(0) = -1, y"Q) = -2e. 
The exact solution for the above problem is y = (1 — x)e x . 

The proposed method is tested on this problem where the domain [0, 1] is divided into 10 equal subintervals. 
The obtained numerical results for this problem are given in Table 1. The maximum absolute error obtained by the 
proposed method is 1.382828xl0~ 5 . 



Table 1: Numerical Results for Example 1 



X 




0.1 


1.192093E-07 


0.2 


3.278255E-06 


0.3 


8.404255E-06 


0.4 


1.186132E-05 


0.5 


1.382828E-05 


0.6 


1.275539E-05 


0.7 


8.046627E-06 


0.8 


3.367662E-06 


0.9 


1.341 105E-07 



Example 2: Consider the linear boundary value problem 

;y (6) + y (5> +sinx y (4> + xy = (2 + sinx + x)e x , 0<;t<l (32) 

subject to 

y(0) = 1, y(l) = e, y'(0) = 1, y'(l) = e, y"(0) = 1, y"(l) = e. 
The exact solution for the above problem is y = e* . 

The proposed method is tested on this problem where the domain [0, 1] is divided into 10 equal subintervals. 
The obtained numerical results for this problem are given in Table 2. The maximum absolute error obtained by the 
proposed method is 1.716614xl0" 5 . 



Table 2: Numerical Results for Example 2 



X 


Absolute Error by the Proposed Method 


0.1 


1.311302E-06 


0.2 


9.536743E-07 


0.3 


4.768372E-07 


0.4 


1.072884E-06 


0.5 


7.510185E-06 
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Table 2: Contd.. 



0.6 


1. 54972 1E-05 


0.7 


1.716614E-05 


0.8 


1.406670E-05 


0.9 


9.775 162E-06 



Example 3: Consider the linear boundary value problem 

y (6) + y'"+y"-y = (-15x 2 +18x-U4)e-\ 0<x<l (33) 
subject to 

y(0) = 0, y(l) = - , y'(0) = 0, y'(l) = - , y"(0) = 0, y"(l) = - . 

e e e 

3 —x 

The exact solution for the above problem is y = X e . 

The proposed method is tested on this problem where the domain [0, 1] is divided into 10 equal subintervals. 
The obtained numerical results for this problem are given in Table 3. The maximum absolute error obtained by the 
proposed method is 2.533 197xl0" 6 . 



Table 3: Numerical Results for Example 3 



X 


Absolute Error by the Proposed Method 


0.1 


3.114110E-08 


0.2 


1.401640E-07 


0.3 


1.527369E-07 


0.4 


3.352761E-08 


0.5 


9.238720E-07 


0.6 


2.115965E-06 


0.7 


2.533197E-06 


0.8 


2.190471E-06 


0.9 


1.639128E-06 



Example 4: Consider the nonlinear boundary value problem 

y (6) + e - x y 2 =e- x + e- 3x , 0<x<l (34) 
subject to 

y(0) = 1, y(l) = -, y'(0) = -1, y'(l) = — , y"(0) = 1, y"(l) = -. 

e e e 

The exact solution for the above problem is y = e x . 

The nonlinear boundary value problem (34) is converted into a sequence of linear boundary value problems 
generated by quasilinearization technique [16] as 

)& + 2e- l y (n) y (n+1) = e-'yf H) + e~ x + e~ 3x , n = 0,1, 2,... (35) 

subject to 
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WO) = °' )Wi>a) = \* W°) = -1, W 1 ) = A W°) = i, W 1 ) =- e 

Here is the (n+l) th approximation for y(x). The domain [0, 1] is divided into 10 equal subintervals and the 

proposed method is applied to the sequence of linear problems (35). The obtained numerical results for this problem are 
presented in Table 4. The maximum absolute error obtained by the proposed method is 2.563000xl0" 6 . 



Table 4: Numerical Results for Example 4 



X 




0.1 


7.748604E-07 


0.2 


4.768372E-07 


0.3 


1.132488E-06 


0.4 


8.940697E-07 


0.5 


1.4901 16E-06 


0.6 


2.563000E-06 


0.7 


2.413988E-06 


0.8 


1.728535E-06 


0.9 


1.102686E-06 



Example 5: Consider the nonlinear boundary value problem 

y (6) +yy (5) -x 3 sm(xx)y'"+yy"+x 2 y 2 = -x 6 cos(xx), 0<x<l (36) 

subject to 

y(0) = 1, y(l) = -1, /(0) = 0, /(l) = 0, /'(0) = -n\ y"{\) = n\ 

The exact solution for the above problem is y = COS(^x). 

The nonlinear boundary value problem (36) is converted into a sequence of linear boundary value problems 
generated by quasilinearization technique [16] as 

(6) t (5) 3 • , x /// // (5) / 

y\ n+ i) + y{n)y\n + i) -x s m (xx)y (n+1) + y {n) y (n+l) + ^ 

+(2x 2 y (n) + y" n) )y(n + i) = y M y"n)+x 2 yl) + y' M yw cos(^) n = 0,1,2,... 

subject to 

J (n+1) (0) = l, y (n+l) (D = -l yU(0) = 0,y[ n+l) (l) = 0, yl +l) {Q) = -x\yl +l) {\) = n 2 . 

Here y( n+ i) is the («+l) th approximation for y(x). The domain [0, 1] is divided into 10 equal subintervals and the 

proposed method is applied to the sequence of linear problems (37). The obtained numerical results for this problem are 
presented in Table 5. The maximum absolute error obtained by the proposed method is 3.288842xl0" 5 . 



Table 5: Numerical Results for Example 5 



X 


Absolute Error by the Proposed Method 


0.1 


4.172325E-07 


0.2 


9.536743E-06 


0.3 


2.193451E-05 


0.4 


3.221631E-05 
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Table 5: Contd., 


0.5 


3.288842E-05 


0.6 


2.413988E-05 


0.7 


1.323223E-05 


0.8 


3.933907E-06 


0.9 


7.748604E-07 



Example 6: Consider the nonlinear boundary value problem 
v (6) -20e- 36v =-40(l + ^r 6 , 0<x<l 



(38) 



subject to 



y(0) = 0, v(l) = y\0) = \, y'(l) = /'(0) = -- y"(l) = ■ 
6 6 12 6 



J_ 
24' 



The exact solution for the above problem is y 



ln(l+jc) 



The nonlinear boundary value problem (38) is converted into a sequence of linear boundary value problems 
generated by quasilinearization technique [16] as 



y (n+1) =120e- 36 ^y (n) + 20e- 3 ^ 



40(1 + x)~ 6 , n = 0,1,2, 



(39) 



subject to 



)W) (°) = °' 0) = » y[n + D (°) = 7 . w = i » >w« (0) = ~ 7 ' y^o d) = -77- 

D D lZ D 24 

Here y(„+;) is the (n+l) lh approximation for y(x). The domain [0, 1] is divided into 10 equal subintervals and the 
proposed method is applied to the sequence of linear problems (39). The obtained numerical results for this problem are 
presented in Table 6. The maximum absolute error obtained by the proposed method is 6.780028xl0" 07 . 

Table 6: Numerical Results for Example 6 



X 


Absolute Error by the Proposed Method Proposed Method 


0.1 


1.303852E-08 


0.2 


1. 67638 1E-08 


0.3 


6.705523E-08 


0.4 


1.192093E-07 


0.5 


3.725290E-07 


0.6 


6.780028E-07 


0.7 


6.854534E-07 


0.8 


5.066395E-07 


0.9 


3.427267E-07 



CONCLUSIONS 

In this paper, we have employed a Petrov-Galerkin method with quintic B-splines as basis functions and septic 
B-splines as weight functions to solve a general sixth order boundary value problems with special case of boundary 
conditions. The quintic B-spline basis set has been redefined into a new set of basis functions which vanish on the 
boundary where the Dirichlet and the Neumann boundary conditions are prescribed. The septic B-splines are redefined into 
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a new set of weight functions which in number match the number of redefined set of basis functions. The proposed method 
has been tested on three linear and three nonlinear sixth order boundary value problems. The numerical results obtained by 
the proposed method are in good agreement with the exact solutions available in the literature. The strength of the 
proposed method lies in its easy applicability, accurate and efficient to solve sixth order boundary value problems. 
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